Source codes: https://github.com/Anne-LiseGerard/dftd_rnaseq

suppressPackageStartupMessages({
    library("reshape2")
    library("gplots")
    library("DESeq2")
    library("mitch")
    library("limma")
    library("kableExtra")
    library("dplyr")
    library("tidyr")
    library("ComplexHeatmap")
    library("RColorBrewer")
    library("circlize")
    library("pathview")
    library("stringr")
    library("grid")
    library("VennDiagram")
    library("enrichplot")
    library("ggplot2")
    library("viridis")
    library("readxl")
    library("cowplot")
    library("grid")
})

knitr::opts_chunk$set(dev = 'svg')

Background

The goal of our study is to compare the expression of genes related to metabolism of DFT1 cell lines to DFT2 cell lines. For this, we have generated transcriptomes for 3 DFT1 and 3 DFT2 cell lines, each in triplicate.

ss <- read.table("../ss.tsv",sep="\t",fill=TRUE,header=TRUE)
ss$DFT <- as.factor(ss$DFT)
ss$clone <- sapply(strsplit(ss$ClientID,"_"),"[[",1)

ss %>%
  kbl(caption="Sample sheet for all samples") %>%
  kable_paper("hover", full_width = F)
Sample sheet for all samples
C ClientID DFT Replicate Index1sequence Index2sequence X.ReadsIdentified.PF. TotalReads TotalBases.Gbases. clone
6180766 4906_1-1 DFT1 1 TTGTTGCA GACGTCGT 2.4465 103481165 15.5 4906
6180767 C5065_1-1 DFT1 1 CCACACTT CGTCATAC 2.8116 118924031 17.8 C5065
6180768 1426_1-1 DFT1 1 CCCGTTTG TTGCTGTA 2.3264 98401219 14.8 1426
6180769 RV_1-1 DFT2 1 ATGCTCCC CCCTGCTG 2.2811 96485136 14.5 RV
6180770 SN_1-1 DFT2 1 GCTCAATA CATTTATC 2.3947 101290147 15.2 SN
6180771 TD549_1-1 DFT2 1 GTAGTTCG CTTGATGC 2.1732 91921221 13.8 TD549
6180772 4906_2-1 DFT1 2 CGAGAACC ATGTATCG 2.2322 94416781 14.2 4906
6180773 C5065_2-1 DFT1 2 GCCATGTA ATAGGGTT 2.4505 103650355 15.5 C5065
6180774 1426_2-1 DFT1 2 TTTCTCTA CTCGACGT 2.4778 104805081 15.7 1426
6180775 RV_2-1 DFT2 2 CCAGCGAT TCTAGTCA 2.9882 126393794 19.0 RV
6180776 SN_2-1 DFT2 2 TGGGAGTG CCCGTCTA 2.4459 103455786 15.5 SN
6180777 TD549_2-1 DFT2 2 CCCTCGTA TAGAAGAG 2.5437 107592495 16.1 TD549
6180778 4906_3-1 DFT1 3 CGATATGG GGGACGTA 2.2711 96062159 14.4 4906
6180779 C5065_3-1 DFT1 3 TTGTGCCC TTTCCATC 2.2763 96282107 14.4 C5065
6180780 1426_3-1 DFT1 3 TGTCCTCT CGACGAAC 2.4952 105541059 15.8 1426
6180781 RV_3-1 DFT2 3 GTATAGTC TTGGTCTC 2.3325 98659234 14.8 RV
6180782 SN_3-1 DFT2 3 TTTGGGAT GTTCACGT 2.6750 113146174 17.0 SN
6180783 TD549_3-1 DFT2 3 CACCAAGC AAAGGGAA 2.1610 91405190 13.7 TD549
DEA5_4NEG DEA4_6NEG NA CGGAGAGG CGAGCGTC 0.0002 8460 0.0 DEA4

We are interested in comparing all DFT1 samples against all DFT2 samples. We will use mSarHar 1.11 from Ensembl v109 for the reference transcriptome.

Functions

# volcanoplots
make_volcano <- function(de,name) {
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$padj),cex=1,pch=19,col="#D5D7E2",
        main=name, xlab="log2 FC", ylab="-log10 pval")
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$padj),cex=1,pch=19,col="#5297A7")
}

# heatmaps
custom_heatmap <- function(zscore, fc, aveexpr, title) {
  h1 <- Heatmap(zscore, cluster_rows = F,
              column_labels = colnames(zscore), name="Z-score",
              cluster_columns = T, 
              column_names_gp = gpar(fontsize = 7, fontface="bold"),
              col=colorRamp2(c(-1.5,0,1.5), hcl_palette = "Blue-Red 2"),
              heatmap_legend_param = list(title = "Z-score", at = seq(-1.5, 1.5,0.5)))
  h2 <- Heatmap(fc, row_labels = rownames(fc), row_names_gp = gpar(fontsize = 6),
              cluster_rows = F, name="logFC", col = col_logFC,column_names_gp = gpar(fontsize =7,fontface="bold"),
              cell_fun = function(j, i, x, y, w, h, col) { # add text to each grid
                grid.text(round(as.numeric(fc[i, j],2),2), x, y,gp=gpar(fontsize=6))
              }, column_title = title, 
              heatmap_legend_param = list(title = "logFC", at = seq(-20, 20, 10)))

  h <- h1 + h2 
  h
}

Load data

Here we load the data in from the aligner.

tmp <- read.table("../fastq/3col.tsv.gz")
x <- as.data.frame(acast(tmp, V2~V1, value.var="V3",fun.aggregate=sum))
dim(x)
## [1] 43512    19

Load gene names.

gn <- read.table("../ref/Sarcophilus_harrisii.mSarHar1.11.cdna+ncrna.genenames.tsv",fill=TRUE)

gn <- gn[order(gn$V1),]

dim(gn)
## [1] 43512     3

Load homology map.

hm <- read.table("../ref/mart_export_ensembl109_2023-07-14.txt",sep="\t",header=TRUE)

Now need to collapse transcript data to genes.

x$gene <- paste(gn$V2,gn$V3)

y <- aggregate(. ~ gene,x,sum)

rownames(y) <- y$gene
y$gene = NULL

dim(y)
## [1] 23829    19

Quality control

Samples with <1M reads should be omitted. Will also round values to integers.

cs <- colSums(y)
cs <- cs[order(cs)]

par(mar=c(5,10,5,2))
barplot(cs,main="All samples",horiz=TRUE,las=1)
abline(v=1e7,col="red",lty=2)

y <- round(y)

MDS

This will help us to visualise the sources of variability in the overall dataset. Fix sample names. Plot MDS of all samples then plot MDS by DFT.

y <- y[,colnames(y) != "DEA5-4NEG"]

ss <- ss[ss$ClientID != "DEA4_6NEG",]

colnames(y) <- sapply(strsplit(ss$ClientID,"-"),"[[",1)

saveRDS(y, file = "y_clines.rds")

cs <- colSums(y)
cs <- cs[order(cs)]

par(mar=c(5,10,5,2))

barplot(cs,main="All samples",horiz=TRUE,las=1, xlab="reads")

cols <- ss$DFT
cols <- gsub("DFT1","#B3B9DF",cols)
cols <- gsub("DFT2","#DDAFB7",cols)
mymds <- plotMDS(y,plot=FALSE)

XMIN=min(mymds$x)
XMAX=max(mymds$x)
XMID=(XMAX+XMIN)/2
XMIN <- XMID + (XMIN-XMID)*1.1
XMAX <- XMID+(XMAX-XMID)*1.1
par(mar = c(5.1, 4.1, 4.1, 2.1) )
plotMDS(mymds,pch=19,cex=3,col=cols,main="Cell lines",xlim=c(XMIN,XMAX))
text(mymds,labels=colnames(y))
mtext("blue=DFT1, pink=DFT2")

y_DFT1 <- y[,colnames(y) %in% c("4906_1","4906_2","4906_3",
                                "1426_1","1426_2","1426_3",
                                "C5065_1","C5065_2","C5065_3")]
y_DFT2 <- y[,colnames(y) %in% c("RV_1","RV_2","RV_3",
                                "SN_1","SN_2","SN_3",
                                "TD549_1","TD549_2","TD549_3")]

par(mar = c(5.1, 4.1, 4.1, 2.1) )
mdsDFT1 <- plotMDS(y_DFT1,plot=FALSE)
plotMDS(mdsDFT1,pch=19,cex=3,col="#B3B9DF",main="MDS plot",xlim=c(XMIN,XMAX))
text(mdsDFT1,labels=colnames(y_DFT1))
mtext("DFT1")

par(mar = c(5.1, 4.1, 4.1, 2.1) )
mdsDFT2 <- plotMDS(y_DFT2,plot=FALSE)
plotMDS(mdsDFT2,pch=19,cex=3,col="#DDAFB7",main="MDS plot",xlim=c(XMIN,XMAX))
text(mdsDFT2, labels=colnames(y_DFT2))
mtext("DFT2")

DESeq2 combining the three replicates

Sum replicates. Plot top 40 DEGs.

x4906 <- rowSums(y[,ss$clone=="4906"])
xC5065 <- rowSums(y[,ss$clone=="C5065"])
x1426 <- rowSums(y[,ss$clone=="1426"])
xRV <- rowSums(y[,ss$clone=="RV"])
xSN <- rowSums(y[,ss$clone=="SN"])
xTD549 <- rowSums(y[,ss$clone=="TD549"])

y2 <- data.frame(x4906,xC5065,x1426,xRV,xSN,xTD549)
ss2 <- as.data.frame(colnames(y2))
colnames(ss2) <- "clone"
ss2$DFT <- factor(c("DFT1","DFT1","DFT1","DFT2","DFT2","DFT2"))

y2 <- y2[which(rowMeans(y2)>10),]
dim(y2) 
## [1] 17468     6
dds <- DESeqDataSetFromMatrix(countData = y2 , colData = ss2, design = ~ DFT )
## converting counts to integer mode
dds2 <- dds
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z<- results(res)
vsd <- vst(dds)
zz<-cbind(as.data.frame(z),assay(vsd))
dge<-as.data.frame(zz[order(zz$pvalue),])
dge2 <- dge
sig <- subset(dge,padj<0.05)
sig2 <- sig
sig2_up <- rownames(subset(sig,log2FoldChange>0))
sig2_dn <- rownames(subset(sig,log2FoldChange<0))
length(sig2_up)
## [1] 4200
length(sig2_dn)
## [1] 4412
# top 20 up & top 20 down
sig <- sig[order(sig$log2FoldChange, decreasing=TRUE),]

sig_noNA <- sig
sig_noNA$genes <- rownames(sig_noNA)
sig_noNA <- filter(sig_noNA, !grepl(' NA', genes))

top <- rbind(head(sig_noNA,20),tail(sig_noNA,20))

mx <- top[,7:ncol(top)] # get normalised counts

mx.scaled <- t(apply(mx[,1:6], 1, scale)) # center and scale each column (Z-score)
colnames(mx.scaled) <- colnames(mx[,1:6])
colnames(mx.scaled) <- c("4906", "C5065", "1426", "RV", "SN", "TD549")
log2FC <- as.matrix(top$log2FoldChange)
rownames(log2FC) <- rownames(top)
colnames(log2FC) <- "logFC"
mean <- as.matrix(top$baseMean)
rownames(mean) <- rownames(top)
colnames(mean) <- "AveExpr"


col_logFC <- colorRamp2(c(-20,0,20), hcl_palette = "Blue-Red 2")

ha <- HeatmapAnnotation(summary=anno_summary(gp=gpar(fill="#D5D7E2"),height=unit(2, "cm")))

rownames(mx.scaled) <- gsub("^.{0,21}", "", rownames(mx.scaled))
rownames(log2FC) <- gsub("^.{0,21}", "", rownames(log2FC))


top40 <- custom_heatmap(mx.scaled, log2FC, mean, title = "Cell lines")

top40

saveRDS(top40, "top40_clines.RDS")
make_volcano(dge2,"Cell lines")

sig[1:50,1:6] %>%
  kbl(caption="Comparison of DFT1 vs DFT2") %>%
  kable_paper("hover", full_width = F)
Comparison of DFT1 vs DFT2
baseMean log2FoldChange lfcSE stat pvalue padj
ENSSHAG00000003963.2 NA 18474.7170 15.79796 1.1976563 13.190727 0.0e+00 0.00e+00
ENSSHAG00000013199.2 PDGFRA 101648.6230 15.67339 0.5200441 30.138574 0.0e+00 0.00e+00
ENSSHAG00000024855.1 SNAI2 3869.1921 15.40132 1.7529564 8.785913 0.0e+00 0.00e+00
ENSSHAG00000004585.2 IBSP 13380.7952 15.33287 1.1634118 13.179224 0.0e+00 0.00e+00
ENSSHAG00000024629.1 NA 7742.8944 14.54383 1.0623179 13.690660 0.0e+00 0.00e+00
ENSSHAG00000030453.1 NA 2021.0210 14.46436 1.2017502 12.036081 0.0e+00 0.00e+00
ENSSHAG00000031194.1 NA 1734.5364 14.24388 1.3557382 10.506367 0.0e+00 0.00e+00
ENSSHAG00000013027.2 RFX4 1314.4713 13.84376 1.9094775 7.250024 0.0e+00 0.00e+00
ENSSHAG00000006714.2 P3H2 2399.3568 13.75004 1.2822571 10.723311 0.0e+00 0.00e+00
ENSSHAG00000004947.2 THBS2 73193.5959 13.50992 0.5132476 26.322420 0.0e+00 0.00e+00
ENSSHAG00000011461.2 PSD4 10786.3556 13.44100 0.6652927 20.203133 0.0e+00 0.00e+00
ENSSHAG00000020999.1 NA 725.7806 12.98694 1.5204866 8.541305 0.0e+00 0.00e+00
ENSSHAG00000028760.1 FGFBP3 9593.7113 12.84541 1.9647464 6.537949 0.0e+00 0.00e+00
ENSSHAG00000015086.2 CYTH4 5920.6132 12.57330 0.6496359 19.354374 0.0e+00 0.00e+00
ENSSHAG00000005458.2 COL3A1 93046.8001 12.38434 0.8065114 15.355438 0.0e+00 0.00e+00
ENSSHAG00000025269.1 NA 465.2029 12.34518 2.2133446 5.577615 0.0e+00 1.00e-07
ENSSHAG00000020986.1 NA 857.2498 12.26509 1.1937914 10.274068 0.0e+00 0.00e+00
ENSSHAG00000003581.2 ATP2B4 3006.5789 12.14268 1.9775155 6.140374 0.0e+00 0.00e+00
ENSSHAG00000002433.2 NA 399.3922 12.12518 1.2167447 9.965262 0.0e+00 0.00e+00
ENSSHAG00000027371.1 NA 1380.2233 12.05697 1.0684268 11.284792 0.0e+00 0.00e+00
ENSSHAG00000015346.2 NA 7887.6771 12.00806 0.7756530 15.481224 0.0e+00 0.00e+00
ENSSHAG00000012585.2 TFAP2B 4877.2129 11.89202 0.8382977 14.185919 0.0e+00 0.00e+00
ENSSHAG00000021326.1 NA 661.5511 11.89100 1.1984707 9.921815 0.0e+00 0.00e+00
ENSSHAG00000005182.2 SERPINE1 21367.7733 11.86676 0.5041184 23.539623 0.0e+00 0.00e+00
ENSSHAG00000008158.2 NA 22953.7196 11.83967 1.5691777 7.545142 0.0e+00 0.00e+00
ENSSHAG00000030465.1 NA 319.5508 11.80335 2.5121349 4.698534 2.6e-06 1.18e-05
ENSSHAG00000003737.2 NA 29315.0122 11.79340 0.3388989 34.799168 0.0e+00 0.00e+00
ENSSHAG00000028636.1 NA 589.0583 11.72350 1.2513273 9.368852 0.0e+00 0.00e+00
ENSSHAG00000006037.2 NA 48768.4674 11.67144 0.7226715 16.150409 0.0e+00 0.00e+00
ENSSHAG00000026258.1 NA 269.9579 11.56028 1.2776044 9.048404 0.0e+00 0.00e+00
ENSSHAG00000017831.2 FAP 523.0322 11.55208 1.6259667 7.104746 0.0e+00 0.00e+00
ENSSHAG00000015915.2 HOXD13 938.6062 11.50197 1.0993823 10.462213 0.0e+00 0.00e+00
ENSSHAG00000009497.2 TNC 28283.5905 11.46547 0.9877163 11.608061 0.0e+00 0.00e+00
ENSSHAG00000024715.1 NA 239.3136 11.38620 1.6966436 6.711013 0.0e+00 0.00e+00
ENSSHAG00000024707.1 NA 460.4332 11.36799 1.2985583 8.754313 0.0e+00 0.00e+00
ENSSHAG00000030987.1 NA 1647.0825 11.31709 0.9252874 12.230887 0.0e+00 0.00e+00
ENSSHAG00000001360.2 NA 432.5737 11.27806 1.1953130 9.435238 0.0e+00 0.00e+00
ENSSHAG00000014465.2 TNR 426.8921 11.25877 1.2275137 9.172014 0.0e+00 0.00e+00
ENSSHAG00000014046.2 CDH15 19589.2177 11.25437 0.5546348 20.291494 0.0e+00 0.00e+00
ENSSHAG00000014110.2 NA 373.3403 11.06549 1.3163357 8.406283 0.0e+00 0.00e+00
ENSSHAG00000025058.1 NA 168.3786 10.87886 1.3334501 8.158432 0.0e+00 0.00e+00
ENSSHAG00000009325.2 BMP5 164.4399 10.84475 1.3301302 8.153148 0.0e+00 0.00e+00
ENSSHAG00000004044.2 ACTRT3 148.5787 10.69836 1.2488369 8.566662 0.0e+00 0.00e+00
ENSSHAG00000022015.1 NA 277.3165 10.63617 1.2310239 8.640099 0.0e+00 0.00e+00
ENSSHAG00000006365.2 CRYBB3 140.9148 10.62197 1.2886846 8.242489 0.0e+00 0.00e+00
ENSSHAG00000002717.2 NA 136.8997 10.58065 1.3696596 7.725020 0.0e+00 0.00e+00
ENSSHAG00000016131.2 GATA2 263.6132 10.56316 1.5752342 6.705769 0.0e+00 0.00e+00
ENSSHAG00000008512.2 CD36 457.5758 10.46304 1.0485317 9.978752 0.0e+00 0.00e+00
ENSSHAG00000016047.2 GALM 447.8834 10.42898 1.1455717 9.103732 0.0e+00 0.00e+00
ENSSHAG00000016786.2 HOXA3 4126.6536 10.38283 0.5189544 20.007206 0.0e+00 0.00e+00
write.table(dge2,file="dge2.tsv",sep="\t")

Enrichment analysis

Need to get human homologs. These were obtained from Ensembl v109 using Biomart.

rownames(dge2) <- sapply(strsplit(rownames(dge2),"\\."),"[[",1)

hm2 <- hm[hm$Tasmanian.devil.gene.stable.ID != "",]

gt <- hm2[,2:3]
length(unique(gt$Tasmanian.devil.gene.stable.ID))
## [1] 16979
length(unique(gt$Gene.name))
## [1] 16646
for(i in 1:length(gt$Gene.name)){
  if(gt$Gene.name[i]==""){gt$Gene.name[i] <- gt$Tasmanian.devil.gene.stable.ID[i]}
}

Now run mitch for DGE. The pathways were sourced from Reactome on 7th July 2023. Plot significant enrichments, then top 20, then pathways within Reactome “metabolism” hierarchy.

genesets <- gmt_import("../ref/ReactomePathways_2023-07-14.gmt")

m2 <- mitch_import(dge2, DEtype="deseq2", geneTable=gt)
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 17468
## Note: no. genes in output = 13563
## Note: estimated proportion of input genes in output = 0.776
head(m2)
##                    x
## A2M     -13.13626596
## A4GALT   -0.07278225
## AAAS     -4.33508400
## AACS      2.13961007
## AADAC    -2.90578338
## AADACL3   1.15402966
res2 <- mitch_calc(m2, genesets, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
if ( !file.exists("mitch_clines_reactome.html") ) {
  mitch_report(res2, "mitch_clines_reactome.html")
}

top <- subset(res2$enrichment_result,p.adjustANOVA<0.05)
saveRDS(top, "~/dftd_RNAseq_annelise/dge/gsea_clines.rds")

top %>%
  kbl(caption="Enriched pathways") %>%
  kable_paper("hover", full_width = F)
Enriched pathways
set setSize pANOVA s.dist p.adjustANOVA
286 Degradation of cysteine and homocysteine 11 0.0004708 -0.6088870 0.0212306
229 Crosslinking of collagen fibrils 16 0.0000992 0.5620709 0.0068156
207 Collagen degradation 47 0.0000000 0.5232821 0.0000004
42 Activation of Matrix Metalloproteinases 22 0.0000217 0.5230344 0.0023319
157 Caspase activation via Death Receptors in the presence of ligand 12 0.0017094 0.5229442 0.0436417
456 GRB2:SOS provides linkage to MAPK signaling for Integrins 12 0.0017300 0.5223600 0.0436417
206 Collagen chain trimerization 30 0.0000011 0.5147073 0.0001397
205 Collagen biosynthesis and modifying enzymes 50 0.0000000 0.4857722 0.0000007
185 Cholesterol biosynthesis 24 0.0000388 0.4852094 0.0032950
171 Cell-extracellular matrix interactions 18 0.0003892 0.4829827 0.0193644
572 Integrin cell surface interactions 63 0.0000000 0.4386596 0.0000006
337 ECM proteoglycans 61 0.0000000 0.4220820 0.0000022
1147 Signal transduction by L1 21 0.0014253 0.4021105 0.0427130
208 Collagen formation 68 0.0000000 0.3860755 0.0000054
573 Integrin signaling 23 0.0014504 0.3836491 0.0427130
1247 Syndecan interactions 24 0.0012870 0.3796440 0.0416276
1243 Sulfur amino acid metabolism 23 0.0021495 -0.3697194 0.0462936
93 Assembly of collagen fibrils and other multimeric structures 47 0.0000156 0.3643719 0.0018740
1235 Smooth Muscle Contraction 29 0.0012237 0.3469831 0.0410639
287 Degradation of the extracellular matrix 107 0.0000000 0.3420029 0.0000005
183 Chemokine receptors bind chemokines 27 0.0021308 0.3415747 0.0462936
1208 Signaling by PDGF 49 0.0001207 0.3176177 0.0075744
1236 Somitogenesis 47 0.0002028 0.3133937 0.0117060
375 Extracellular matrix organization 231 0.0000000 0.3023319 0.0000000
776 Non-integrin membrane-ECM interactions 49 0.0004355 0.2906132 0.0202699
872 Platelet degranulation 88 0.0000259 0.2596770 0.0024947
1019 Regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins (IGFBPs) 85 0.0000365 0.2592949 0.0032927
1074 Response to elevated platelet cytosolic Ca2+ 93 0.0000697 0.2388941 0.0052913
1133 Semaphorin interactions 56 0.0020537 0.2382812 0.0462936
418 Formation of paraxial mesoderm 60 0.0020070 0.2307463 0.0459706
883 Post-translational protein phosphorylation 76 0.0007652 0.2234588 0.0306723
347 ER-Phagosome pathway 77 0.0007355 0.2227312 0.0303231
368 Eukaryotic Translation Elongation 67 0.0016960 0.2219276 0.0436417
546 Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell 71 0.0012982 0.2209405 0.0416276
414 Formation of a pool of free 40S subunits 76 0.0011072 0.2166004 0.0395156
77 Antigen processing-Cross presentation 87 0.0006988 0.2104900 0.0296577
648 MAPK6/MAPK4 signaling 78 0.0013520 0.2100834 0.0423523
595 Interleukin-4 and Interleukin-13 signaling 78 0.0016034 0.2068443 0.0436417
458 GTP hydrolysis and joining of the 60S ribosomal subunit 87 0.0010061 0.2042072 0.0382042
619 L13a-mediated translational silencing of Ceruloplasmin expression 86 0.0011228 0.2034474 0.0395156
311 Diseases of glycosylation 106 0.0008104 0.1885306 0.0316073
149 Cap-dependent Translation Initiation 93 0.0017700 0.1877977 0.0436417
369 Eukaryotic Translation Initiation 93 0.0017700 0.1877977 0.0436417
771 Neutrophil degranulation 351 0.0000000 0.1718105 0.0000054
870 Platelet activation, signaling and aggregation 193 0.0000599 0.1679516 0.0048003
1215 Signaling by ROBO receptors 172 0.0019622 0.1371500 0.0456688
445 G2/M Transition 175 0.0018011 0.1370959 0.0436417
1359 Transmission across Chemical Synapses 195 0.0010384 -0.1365791 0.0384190
708 Mitotic G2-G2/M phases 177 0.0018428 0.1360347 0.0436417
709 Mitotic Metaphase and Anaphase 206 0.0016564 0.1275027 0.0436417
565 Innate Immune System 718 0.0000000 0.1267316 0.0000021
706 Mitotic Anaphase 205 0.0021115 0.1248964 0.0462936
927 RHO GTPase Effectors 232 0.0012091 0.1237380 0.0410639
524 Hemostasis 447 0.0000226 0.1176378 0.0023319
640 M Phase 326 0.0003304 0.1161895 0.0176607
1347 Translation 254 0.0015502 0.1157384 0.0436417
1188 Signaling by Interleukins 343 0.0013795 0.1009986 0.0423523
545 Immune System 1394 0.0000000 0.0978397 0.0000006
103 Axon guidance 467 0.0004197 0.0958918 0.0201878
1218 Signaling by Rho GTPases 570 0.0001044 0.0958489 0.0068487
59 Adaptive Immune System 585 0.0000951 0.0952178 0.0068156
1219 Signaling by Rho GTPases, Miro GTPases and RHOBTB3 585 0.0001616 0.0920456 0.0097133
764 Nervous system development 492 0.0018449 0.0825619 0.0436417
236 Cytokine Signaling in Immune system 514 0.0017094 0.0814258 0.0436417
296 Developmental Biology 790 0.0004940 0.0737359 0.0215998
300 Disease 1367 0.0003517 0.0588440 0.0181266
1144 Signal Transduction 1971 0.0002946 0.0509136 0.0163530
ggplot(top, aes(x = s.dist, y = reorder(set, s.dist), color = p.adjustANOVA, size = setSize)) + 
  geom_point(stat = 'identity') + 
  scale_color_viridis(option="mako", direction=1,limits = c(0, 0.05)) +
  xlab("s dist") + ylab("pathway") + ggtitle("Reactome ALL") + 
  theme_bw()

top20 <- rbind(slice_max(top, order_by=s.dist, n=10), slice_min(top, order_by=s.dist, n=10))

ggplot(top20, aes(x = s.dist, y = reorder(set, s.dist), color = p.adjustANOVA, size = setSize)) + 
  geom_point(stat = 'identity') + 
  scale_color_viridis(option="mako", direction=1,limits = c(0, 0.05)) +
  xlab("s dist") + ylab("pathway") + ggtitle("Reactome top 20") + 
  theme_bw()

react_metabosets <- readRDS("../ref/reactome_metabo.rds")

metabo1 <- filter(top, set %in% names(react_metabosets))
saveRDS(metabo1, "~/dftd_RNAseq_annelise/dge/gsea_clines_metabo.rds")

ggplot(metabo1 , aes(x = s.dist, y = reorder(set, s.dist), color = p.adjustANOVA, size = setSize)) + 
  geom_point(stat = 'identity') + 
  scale_color_viridis(option="mako", direction=1,limits = c(0, 0.05)) +
  xlab("s dist") + ylab("pathway") + ggtitle("Reactome metabolism") + 
  theme_bw()

top_reactome <- top

Visualisation

Heatmaps. First plot enriched pathways, then plot pathways of interest for which experimental data is available.

# see https://www.youtube.com/watch?v=ht1r34-ifVI for reference

hm3 <- data.frame(ENSSHAG=hm2[,2],geneID=hm2[,3])
dge2 <- dge
dge2$ensembl <- rownames(dge2)
split <- strsplit(dge2$ensembl, split = " ")
ensembl <- c()
gene_name <- c()
for (i in 1:length(split)){
  ensembl <- append(ensembl, split[[i]][1]) 
  gene_name <- append(gene_name, split[[i]][2]) 
}
dge2$ensembl <- ensembl
dge2$gene_name <- gene_name

dge3 <- dge2

# only plot significant DEGs
dge3sig <- filter(dge3, padj < 0.05)
dge3sig <- dge3sig[order(dge3sig$log2FoldChange, decreasing=TRUE),]

mx <- dge3sig[,7:ncol(dge3sig)] # get normalised counts
 
mx.scaled <- t(apply(mx[,1:6], 1, scale)) # center and scale each column (Z-score)
colnames(mx.scaled) <- c("DFT1 4906", "DFT1 C5065", "DFT1 1426", "DFT2 RV", "DFT2 SN", "DFT2 TD549")
log2FC <- as.matrix(dge3sig$log2FoldChange)
rownames(log2FC) <- rownames(dge3sig)
colnames(log2FC) <- "logFC"
mean <- as.matrix(dge3sig$baseMean)
rownames(mean) <- rownames(dge3sig)
colnames(mean) <- "AveExpr"

col_logFC <- colorRamp2(c(-20,0,20), hcl_palette = "Blue-Red 2")

custom_heatmap <- function(zscore, fc, aveexpr, title) {
  h1 <- Heatmap(zscore, cluster_rows = F,
              column_labels = colnames(zscore), name="Z-score",
              cluster_columns = T, 
              column_names_gp = gpar(fontsize = 7, fontface="bold"),
             col=colorRamp2(c(-1.5,0,1.5), hcl_palette = "Blue-Red 2"),
             heatmap_legend_param = list(title = "Z-score", at = seq(-1.5, 1.5, 0.5)))
  h2 <- Heatmap(fc, row_labels = rownames(fc), row_names_gp = gpar(fontsize = 6),
              cluster_rows = F, name="logFC", col = col_logFC,column_names_gp = gpar(fontsize =7,fontface="bold"),
              cell_fun = function(j, i, x, y, w, h, col) { # add text to each grid
                grid.text(round(as.numeric(fc[i, j],2),2), x, y,gp=gpar(fontsize=6))
              }, column_title = title, 
              heatmap_legend_param = list(title = "logFC", at = seq(-20, 20, 10)))
  
  h <- h1 + h2 
  h
}

ha <- HeatmapAnnotation(summary=anno_summary(gp=gpar(fill="#D5D7E2"),height=unit(1.5, "cm")))

# ----- Top DEGs in REACTOME metabolic pathways ------

metabo <- c() # get all names of genes involved in metabolism
for(i in 1:length(genesets)){metabo <- c(metabo,genesets[[i]])}
metabo <- unique(metabo) # get rid of duplicates

temp <- c() # get ENSSHAG and gene ID
for(i in 1:length(metabo)){temp <- rbind(temp,filter(hm3, geneID==metabo[i]))}
metabolism <- c(t(unite(temp,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

# collapse transcript IDs to gene IDs
rownames(mx.scaled) <- sub('\\..*\\s', ' ', rownames(mx.scaled))
rownames(log2FC) <- sub('\\..*\\s', ' ', rownames(log2FC))
rownames(mean) <- sub('\\..*\\s', ' ', rownames(mean))
dge3$ensembl <- sub('\\..*', '', dge3$ensembl)
rownames(dge3) <- sub('\\..*\\s', ' ', rownames(dge3))

rownames(mx.scaled) <- gsub("^.{0,19}", "", rownames(mx.scaled))
rownames(log2FC) <- gsub("^.{0,19}", "", rownames(log2FC))

metabolism <- gsub("^.{0,19}", "", metabolism)

met <- subset(mx.scaled, rownames(mx.scaled) %in% metabolism)
metFC <- subset(log2FC, rownames(log2FC) %in% metabolism)
metAE <- subset(mean, rownames(mean) %in% metabolism)

# top 20 up & down
met40 <- rbind(head(met,20),tail(met,20))
met40FC <- as.matrix(c(head(metFC,20),tail(metFC,20)))
rownames(met40FC) <- rownames(met40)
colnames(met40FC) <- "logFC"
met40AE <- as.matrix(c(head(metAE,20),tail(metAE,20)))

custom_heatmap(met40, met40FC, met40AE, title = "TOP 40 - REACTOME metabolism")

# ----- All DEGs in REACTOME metabolic pathways ------

dge3_metabo_reactome <- filter(dge3, gene_name %in% metabo)

make_volcano(dge3_metabo_reactome, "Cell lines - Metabolism genes - REACTOME")

# ----- Individual TOP REACTOME metabolic pathways -----

pathways <- metabo1$set
heatmaps <- vector(mode = "list", length = length(pathways))
for(i in 1:length(pathways)){
  
  metabo <- c() # get all names of genes involved in metabolism
  for(j in 1:length(genesets)){metabo <- c(metabo,genesets[[pathways[i]]])}
  metabo <- unique(metabo) # get rid of duplicates

  temp <- c() # get ENSSHAG and gene ID
  for(k in 1:length(metabo)){temp <- rbind(temp,filter(hm3, geneID==metabo[k]))}
  metabolism <- c(t(unite(temp,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))
  
  metabolism <- gsub("^.{0,19}", "", metabolism)
  
  met <- subset(mx.scaled, rownames(mx.scaled) %in% metabolism)
  metFC <- subset(log2FC, rownames(log2FC) %in% metabolism)

  heatmaps[[i]] <- custom_heatmap(met, metFC, title = paste(pathways[i]))
  
}
print(heatmaps[[1]])

saveRDS(heatmaps[[1]], "cyst.RDS")
print(heatmaps[[2]])

saveRDS(heatmaps[[2]], "chol.RDS")
print(heatmaps[[3]])

saveRDS(heatmaps[[3]], "sulf.RDS")
# ----- Other interesting pathways REACTOME -----

mm <- filter(hm3, geneID %in% genesets[["Glycolysis"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)
metAE <- subset(mean, rownames(mean) %in% mm)

glyco_p <- custom_heatmap(met, metFC, title = "Cell lines")

glyco_p

saveRDS(glyco_p, "glyco_p.rds")
mm <- filter(hm3, geneID %in% genesets[["Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins."]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)
metAE <- subset(mean, rownames(mean) %in% mm)

custom_heatmap(met, metFC, metAE, title = "Respiratory electron transport, \n ATP synthesis by chemiosmotic coupling, \n and heat production by uncoupling proteins.")

mm <- filter(hm3, geneID %in% genesets[["Citric acid cycle (TCA cycle)"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)
metAE <- subset(mean, rownames(mean) %in% mm)

custom_heatmap(met, metFC, metAE, title = "Citric acid cycle (TCA cycle)")

mm <- filter(hm3, geneID %in% genesets[["Respiratory electron transport"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)
metAE <- subset(mean, rownames(mean) %in% mm)

custom_heatmap(met, metFC, metAE, title = "Respiratory electron transport")

mm <- filter(hm3, geneID %in% genesets[["DNA Repair"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

custom_heatmap(met, metFC, title = "DNA Repair")

mm <- filter(hm3, geneID %in% genesets[["DNA Double-Strand Break Repair"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

custom_heatmap(met, metFC, title = "DNA Double-Strand Break Repair")

mm <- filter(hm3, geneID %in% genesets[["Detoxification of Reactive Oxygen Species"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

ROS_p <- custom_heatmap(met, metFC, title = "Cell lines")

ROS_p 

saveRDS(ROS_p, "ROS_p.rds")
mm <- filter(hm3, geneID %in% genesets[["Metabolism of nucleotides"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

nucl_p <- custom_heatmap(met, metFC, title = "Metabolism of nucleotides")

nucl_p 

saveRDS(nucl_p, "nucl.rds")
mm <- filter(hm3, geneID %in% genesets[["Pentose phosphate pathway"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

custom_heatmap(met, metFC, title = "Pentose phosphate pathway")

mm <- filter(hm3, geneID %in% genesets[["Fatty acid metabolism"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

custom_heatmap(met, metFC, title = "Fatty acid metabolism")

mm <- filter(hm3, geneID %in% genesets[["Glutamate and glutamine metabolism"]])
mm <- c(t(unite(mm,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

mm <- gsub("^.{0,19}", "", mm)

met <- subset(mx.scaled, rownames(mx.scaled) %in% mm)
metFC <- subset(log2FC, rownames(log2FC) %in% mm)

custom_heatmap(met, metFC, metAE, title = "Glutamate and glutamine metabolism")

Plotting custom gene sets

Two custom gene sets were added to the previous analysis containing cell competition & metastasis related genes based on a literature search.

cell_compet_genes <- read_excel("../ref/gsea/Cell_compet_genes.xlsx")
names(cell_compet_genes)[names(cell_compet_genes) == "Ensembl_ID"] <- "Gene.stable.ID"

cell_compet_genes <- inner_join(cell_compet_genes,hm,by="Gene.stable.ID") # making sure these genes exist in devils
## Warning in inner_join(cell_compet_genes, hm, by = "Gene.stable.ID"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 27 of `x` matches multiple rows in `y`.
## ℹ Row 65234 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
genesets$"Cell Competition" <- unique(cell_compet_genes$Gene.name)
genesets$"Cell Competition - Mechanical Competition" <- unique(filter(cell_compet_genes, Type == "Mechanical_Competition")$Gene.name)

metastasis_genes <- read_excel("../ref/gsea/Metastasis_genes.xlsx")
names(metastasis_genes)[names(metastasis_genes) == "Ensembl_ID"] <- "Gene.stable.ID"

metastasis_genes <- inner_join(metastasis_genes,hm,by="Gene.stable.ID") # making sure these genes exist in devils

genesets$"Metastasis" <- unique(metastasis_genes$Gene.name)
genesets$"Metastasis - promotors" <- unique(filter(metastasis_genes, Type == "Metastasis_promotor")$Gene.name)
genesets$"Metastasis - suppressors" <- unique(filter(metastasis_genes, Type == "Metastasis_supressor")$Gene.name)

saveRDS(genesets, "~/dftd_RNAseq_annelise/ref/gsea/reactome_plus_custom.RDS")

res_custom <- mitch_calc(m2, genesets, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
if ( !file.exists("mitch_custom.html") ) {
  mitch_report(res_custom, "mitch_custom.html")
}

cc <- filter(hm3, geneID %in% genesets[["Cell Competition"]])
cc <- c(t(unite(cc,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

cc <- gsub("^.{0,19}", "", cc)

met <- subset(mx.scaled, rownames(mx.scaled) %in% cc)
metFC <- subset(log2FC, rownames(log2FC) %in% cc)

custom_heatmap(met, metFC, title = "Cell competition")

meta <- filter(hm3, geneID %in% genesets[["Metastasis"]])
meta <- c(t(unite(meta,"ENSSHAG_geneID", ENSSHAG:geneID, sep=" ")))

meta <- gsub("^.{0,19}", "", meta)

met <- subset(mx.scaled, rownames(mx.scaled) %in% meta)
metFC <- subset(log2FC, rownames(log2FC) %in% meta)

custom_heatmap(met, metFC, metAE, title = "Metastasis")

Session information

For reproducibility.

sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
##  [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
##  [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.3               readxl_1.4.3               
##  [3] viridis_0.6.5               viridisLite_0.4.2          
##  [5] ggplot2_3.5.1               enrichplot_1.24.4          
##  [7] VennDiagram_1.7.3           futile.logger_1.4.3        
##  [9] stringr_1.5.1               pathview_1.44.0            
## [11] circlize_0.4.16             RColorBrewer_1.1-3         
## [13] ComplexHeatmap_2.20.0       tidyr_1.3.1                
## [15] dplyr_1.1.4                 kableExtra_1.4.0           
## [17] limma_3.60.4                mitch_1.16.0               
## [19] DESeq2_1.44.0               SummarizedExperiment_1.34.0
## [21] Biobase_2.64.0              MatrixGenerics_1.16.0      
## [23] matrixStats_1.4.1           GenomicRanges_1.56.1       
## [25] GenomeInfoDb_1.40.1         IRanges_2.38.1             
## [27] S4Vectors_0.42.1            BiocGenerics_0.50.0        
## [29] gplots_3.1.3.1              reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.4.2           later_1.3.2             ggplotify_0.1.2        
##   [4] bitops_1.0-8            cellranger_1.1.0        tibble_3.2.1           
##   [7] R.oo_1.26.0             polyclip_1.10-7         graph_1.82.0           
##  [10] XML_3.99-0.17           lifecycle_1.0.4         httr2_1.0.3            
##  [13] doParallel_1.0.17       lattice_0.22-6          MASS_7.3-61            
##  [16] magrittr_2.0.3          sass_0.4.9              rmarkdown_2.28         
##  [19] jquerylib_0.1.4         yaml_2.3.10             httpuv_1.6.15          
##  [22] DBI_1.2.3               abind_1.4-5             zlibbioc_1.50.0        
##  [25] purrr_1.0.2             R.utils_2.12.3          ggraph_2.2.1           
##  [28] RCurl_1.98-1.16         yulab.utils_0.1.7       tweenr_2.0.3           
##  [31] rappdirs_0.3.3          GenomeInfoDbData_1.2.12 ggrepel_0.9.6          
##  [34] tidytree_0.4.6          svglite_2.1.3           codetools_0.2-20       
##  [37] DelayedArray_0.30.1     DOSE_3.30.5             xml2_1.3.6             
##  [40] ggforce_0.4.2           tidyselect_1.2.1        shape_1.4.6.1          
##  [43] aplot_0.2.3             UCSC.utils_1.0.0        farver_2.1.2           
##  [46] jsonlite_1.8.8          GetoptLong_1.0.5        tidygraph_1.3.1        
##  [49] iterators_1.0.14        systemfonts_1.1.0       foreach_1.5.2          
##  [52] tools_4.4.2             treeio_1.28.0           Rcpp_1.0.13            
##  [55] glue_1.7.0              gridExtra_2.3           SparseArray_1.4.8      
##  [58] xfun_0.49               qvalue_2.36.0           withr_3.0.1            
##  [61] formatR_1.14            fastmap_1.2.0           GGally_2.2.1           
##  [64] fansi_1.0.6             caTools_1.18.3          digest_0.6.37          
##  [67] gridGraphics_0.5-1      R6_2.5.1                mime_0.12              
##  [70] colorspace_2.1-1        Cairo_1.6-2             GO.db_3.19.1           
##  [73] gtools_3.9.5            RSQLite_2.3.7           R.methodsS3_1.8.2      
##  [76] utf8_1.2.4              generics_0.1.3          data.table_1.16.0      
##  [79] graphlayouts_1.1.1      httr_1.4.7              htmlwidgets_1.6.4      
##  [82] S4Arrays_1.4.1          scatterpie_0.2.4        ggstats_0.6.0          
##  [85] pkgconfig_2.0.3         gtable_0.3.5            blob_1.2.4             
##  [88] XVector_0.44.0          shadowtext_0.1.4        htmltools_0.5.8.1      
##  [91] fgsea_1.30.0            echarts4r_0.4.5         clue_0.3-65            
##  [94] scales_1.3.0            png_0.1-8               ggfun_0.1.6            
##  [97] knitr_1.48              lambda.r_1.2.4          rstudioapi_0.16.0      
## [100] rjson_0.2.22            nlme_3.1-166            org.Hs.eg.db_3.19.1    
## [103] cachem_1.1.0            GlobalOptions_0.1.2     KernSmooth_2.23-24     
## [106] parallel_4.4.2          AnnotationDbi_1.66.0    pillar_1.9.0           
## [109] vctrs_0.6.5             promises_1.3.0          xtable_1.8-4           
## [112] cluster_2.1.6           beeswarm_0.4.0          Rgraphviz_2.48.0       
## [115] evaluate_0.24.0         KEGGgraph_1.64.0        magick_2.8.4           
## [118] cli_3.6.3               locfit_1.5-9.10         compiler_4.4.2         
## [121] futile.options_1.0.1    rlang_1.1.4             crayon_1.5.3           
## [124] labeling_0.4.3          plyr_1.8.9              fs_1.6.4               
## [127] stringi_1.8.4           BiocParallel_1.38.0     munsell_0.5.1          
## [130] Biostrings_2.72.1       lazyeval_0.2.2          GOSemSim_2.30.2        
## [133] Matrix_1.7-0            patchwork_1.2.0         bit64_4.0.5            
## [136] KEGGREST_1.44.1         statmod_1.5.0           shiny_1.9.1            
## [139] highr_0.11              igraph_2.0.3            memoise_2.0.1          
## [142] bslib_0.8.0             ggtree_3.12.0           fastmatch_1.1-4        
## [145] bit_4.0.5               ape_5.8